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The 4-th order Runge-Kutta method in the complex plane is proposed for numerically advancing the solutions 
of a system of first order differential equations in one external invariant satisfied by the master integrals related to 
a Feynman graph. The particular case of the general massive 2-loop sunrise self-mass diagram is analyzed. The 
method offers a reliable and robust approach to the direct and precise numerical evaluation of master integrals. 


1. Introduction 

The very high precision of present and planned 
particle physics experiments requires comparable 
or better accuracy on the theoretical side. This 
fact promotes developments of new methods in 
the calculations of radiative corrections, which 
are today a living and expanding field. 

The nowadays widespread organization of the 
calculations is based on the integration by part 
identities and on the evaluation of the master in¬ 
tegrals (MI) [1]. We believe that the systematic 
use of the differential equations for the MI, or 
Master Differential Equations (MDE), can be a 
viable method for their analytic calculations in 
many cases. In these cases, but also when the 
number of variables and parameters prevents the 
success of an analytic calculation, the MDE can 
still be profitably used for direct numerical eval¬ 
uation of the MI. This is an alternative to the 
more commonly used integration methods or to 
the more recently introduced difference equations 
method. 

A method which uses the MDE to get a nu¬ 
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merical solution, starting from a known value, is 
presented here and its features are discussed. 


2. Master Differential Equations 


Starting from the integral representation of the 
Nmi MI, related to a certain Feynman graph, 
by derivation with respect to one of the internal 
masses rrii [2] or one of the external invariants s e 
[3] and with the repeated use of the integration 
by part identities, a system of Nmi independent 
first order partial MDE is obtained for the Nmi 
MI. For any of the s e , say Sj, the equations have 
in general the form 


K k {r 




.) = 


X M k,i(n, m?,s e )Fi(n, m 2 , s e ) + T k {n, m 2 , s e ), 

i 


k,l = 1,..., Nmi 


(1) 


where F k (n,m 2 ,s e ) are the MI, K k (m 2 ,s e ) 
and M k: i(n,m?,s e ) are polynomials, while 
T k (n,m 2 , s e ) are polynomials times simpler MI 
of the subgraphs of the considered graph. The 
roots of the equations 

K k (m 2 ,s e ) = 0 (2) 


identify the special points, where numerical cal¬ 
culations are troublesome. Fortunately analytic 
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calculations at those points come out to be possi¬ 
ble in all the attempted cases so far. They might 
not be simple and often require some external 
knowledge, like the assumption of regularity of 
the solution at that special point. 

To solve the system of equations it is neces¬ 
sary to know the MI for a chosen value of the 
differential variable, Sj in Eq.(l). For that pur¬ 
pose we use the analytic expressions at the spe¬ 
cial points, taken as the starting points of the 
advancing solution path. Moreover starting from 
one special point, not only the values of the MI 
are necessary, but also their first order deriva¬ 
tives at that point. That is because some of the 
coefficients Kk(m 2 ,s e ) of the MI derivatives in 
the differential equations Eq.(l) vanish at that 
point. Therefore also the analytic expressions for 
the first derivatives of MI at special points are ob¬ 
tained, but this usually comes out to be a simpler 
task (unless poles in the limit of the number of 
dimensions n going to 4 are present). 

Enlarging the number of loops and legs in¬ 
creases the number of parameters, MI and equa¬ 
tions, but does not change or spoil the method. 


3. The 4-th order Runge-Kutta method 


Many methods are available for obtaining the 
numerical solutions of a first-order differential 
equation [4] 


dy{x) 

dx 


f(x,y) . 


( 3 ) 


The Euler method advances the solution from a 
point x n . where the solution y n is known, to the 
point x n+ \ = x n + A 


Vn+I =Vn + Af(x n ,y n ) + 0(A 2 ) (4) 

omitting terms of order A 2 . A direct improve¬ 
ment of the Euler method is the 4-th order Runge- 
Kutta method, that we choose, because it is con¬ 
sidered a rather precise and robust approach. By 
suitably choosing the intermediate points where 
calculating f(x,y) one obtains the 4-th order 
Runge-Kutta formula 


ki = A f(x n ,y n ), 

i . A k u 

k 2 = A }{x n + —,y n + —), 


h = Af(x n + y,2M + y), 

h = A/(a : n + A, y n + k 3 ), 

ki k 2 k 3 /c4 . c. , . 

Vn+l = 3/tt. H ^—I g—I ^—I—^—b 0 {A ) (5) 

which omits terms of order A 5 . 

To avoid numerical problems due to the pres¬ 
ence of special points on the real axis, it is conve¬ 
nient to choose a path for advancing the solution 
in the complex plane of x. 

The extension from one first-order differential 
equation to a system of Nmi first-order MDE for 
the Nmi MI is straightforward [4], 

4. Results: sunrise, ... 

To test the method we have chosen to start 
from the simple, but not trivial, 2-loop sunrise 
graph with arbitrary masses [5,6], shown in Fig.l. 



Figure 1. The general massive 2-loop sunrise self¬ 
mass diagram. 


This graph is one of the topologies of the 2-loop 
self-mass and has 4 MI. The other topologies with 
4 and 5 propagators, shown in Fig.2 and in Fig.3 
respectively, have one more MI each [7,5,8]. The 
calculations for these graphs are in progress [9]. 

The sunrise general amplitude can be written 
in the integral form as 


A(n,m 2 ,p 2 , -ai, -a 2 ,-a 3 , fit,/3 2 ) = 


r d n k 1 r d n k 2 (p. fci)ft(p- k 2 )fo 
J (27r)"-2 J (27r)™- 2 D* 1 L>“ 2 Z?“ 3 


D\ — (ki + to 2 ), 

D 2 = (fea+TOj), 
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Figure 2. The general massive 2-loop 4- 
denominators self-mass diagram. 



Figure 3. The general massive 2-loop 5- 
denominators self-mass diagram. 

D 3 = {{p- h- k 2 ) 2 + m§), (6) 

where p is the arbitrary mass scale accounting for 
the continuous value of the dimensions n. As the 
natural scale of the problem is the threshold of the 
sunrise amplitudes, we take p = m\ + m 2 + m 3 . 
We choose for the 4 MI the amplitudes 

F 0 (n,m 2 ,p 2 ) =A(n,mf,p 2 ,- 1,-1,-1,0,0) , 
Fi(n,m 2 ,p 2 ) =A(n,m 2 ,p 2 ,- 2,-1,-1,0,0) , 
F 2 (n,m 2 ,p 2 ) =A(n,m 2 ,p 2 ,- 1,-2,-1,0,0) , 
F 3 (n,m 2 ,p 2 ) =A(n,m 2 ,p 2 ,- 1,-1,-2,0,0) , (7) 

which are connected by the relations (j, i = 
1,2,3) 

Fj(n, m 2 ,p 2 ) = F o(n,m 2 ,p 2 ) . (8) 

The above relations are very important for ana¬ 
lytic calculations, but are not used in the present 
numerical approach. 

To obtain the system of differential equations 
for the 4 MI in the form of Eq.(l) it is neces¬ 
sary to develop explicitly the integration by parts 


identities for the amplitudes with the values of 
the exponents ctj and f3j satisfying the relations 
1 , 2,3 ( a * — 1) = 2 and Y2j= 1,2 = 2- 

In the differential equations of the sunrise MI 
the only lower order diagram entering is the 1- 
loop vacuum graph 

rru 2\ /' d n k 1 

lin.m ) = / 7 —^— 7—7 - 7 

J (27r) ra_2 k 2 + m 2 

m n ~ 2 C(n) 

(n — 2 )(n — 4) 

The function 

C(n) = (2^) (4 ' n) r(3-|) , 

which appears in the expressions for the MI as 
an overall factor with an exponent equal to the 
number of loops, is usually kept unexpanded in 
the limit n —■> 4, and only at the very end of the 
calculation for finite quantities is set (7(4) = 1. 

When the sunrise MI are expanded in (?i — 4), 
for j = 0,1, 2, 3, and * = 1,2, 3, 


(9) 

( 10 ) 


Fj{n,m 2 ,p 2 ) = C 2 (n) 


1 


1 


Fj 2) (™?,P 2 ) 


(n — 4 ) 2 0 


(n _ 4 fj 1 } K’P 2 ) + F j°\m 2 , p 2 ) 


+0(n — 4) 




( 11 ) 


the coefficients of the poles can be easily obtained 
analytically for arbitrary values of the external 
squared momentum p 2 , 

F o~ 2 \m 2 ,p 2 ) = -\{m\ + ml + m§) , 


F 


— 1 j P 2 


’(m 2 ,p 2 ) = + ^(ml+m 2 2 + ml) 




F i 2 \rn 2 ,p 2 ) = - , k = 1,2,3 

= ~ + 1 Mf) ■ 

The finite parts satisfy the differential equations 
P 2 ^F'o ) {m 2 i ,p 2 )= F ^\m 2 ,p 2 ) + -Po _ 1 ) (m?,p 2 ) 
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+ X! m2 j F j°\ m i’P 2 ) , (!3) 

J=l,2,3 


and ( i,j, k, l = 1, 2, 3, with j ^ k l) 


8D(m 2 , p 2 )p 2 -^~oF^ 0) (m 2 , p 2 ) = 


d_ 

dp 2 

4 ■D(m 2 i ,p 2 )F < /~ 1) (m 2 ,p 2 ) 

+ Pip(m 2 ,p 2 ) [l6F 0 (0) (m 2 ,p 2 ) + 28*’ 0 ( - 1) (m?,p 2 ) 

+12F 0 ( - 2) (m 2 ,p 2 ) 

+ 8 Pi,i(m 2 ,p 2 ) F[ 0) (m 2 ,p 2 ) + f[~ 1] (m 2 ,p 2 ) 

+ 8 Pi d (mlp 2 ) [F^(m 2 ,p 2 ) + ^(mf.p 2 ) 

+ 8P/, fe (m 2 ,p 2 ) F^ 0 \m 2 ,p 2 ) + P^ _1) (m 2 ,p 2 ) 


+ Qi,i{m 2 ,p 2 ) m 2 m 2 [log(m 2 ) + log(m £)] 2 
+ Qi,j(m 2 ,p 2 ) mfm 2 k [log(rn 2 ) + log(m 2 )] 2 
+ Qi,k(m 2 ,p 2 ) ?n^m 2 [log(m 2 ) + log(m 2 )] 2 . (14) 
The special points are p 2 = 0, oo and the roots of 


D{m 2 ,p 2 ) = [p 2 + (mi + ?n 2 + m 3 ) 2 ] 

[p 2 + (mi + m 2 - m 3 ) 2 ] [p 2 + (mi - m 2 + m 3 ) 2 ] 
[p 2 + (mi - m 2 - m 3 ) 2 ] = 0 , (15) 

and Pij(m 2 ,p 2 ) and Qi,j(m^,p 2 ) are polynomials 
in p 2 and in the masses, whose explicit expres¬ 
sions can be found in [5]. 

From these equations the analytic expressions 
for their first order expansion were completed 
around the special points [5,10,11,6]: p 2 = 0; 
p 2 = oo; p 2 = —(mi + m 2 + m3) 2 , the threshold; 
p 2 = —(mi + m 2 — m3) 2 , the pseudo-thresholds. 

To obtain numerical results for arbitrary values 
of p 2 , a 4th-order Runge-Kutta formula is imple¬ 
mented in a FORTRAN code, with a solution ad¬ 
vancing path starting from the special points, so 
that also the first term in the expansion is neces¬ 
sary. 

The path followed starts usually from p 2 = 0 
and moves in the lower half complex plane of 
p 2 = p 2 /(m\ + m 2 + m3) 2 , as shown in Fig.4, to 
avoid proximity to the other special points, which 
can cause loss in precision. Values between special 
points can be safely reached through a complex 


path as also shown in Fig.4. For values ofp 2 very 
close to a special point, we start from the analyt¬ 
ical expansion at that special point. Subtracted 
differential equations are used when starting from 
p 2 = 00 or from threshold, as that points are not 
regular points of the MDE Eq.(13),Eq.(14). 

Remarkable self-consistency checks are easily 
provided by comparing the results obtained ei¬ 
ther starting from the same point and choosing 
different paths to arrive to the same final point, 
or choosing directly different starting points and 
again arriving to the same final point. 
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Imp^ 


-1 
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+ 


Re p r 2 


- 0.1 


Figure 4. Paths followed in the complex p 2 plane. 
On the real axis are indicated the positions of the 
threshold (- 1 ) and pseudo-thresholds. 


The execution of the program is rather fast and 
precise: with an Intel Pentium III of 1 GHz we 
get values with 7 digits requiring times ranging 
from a fraction of a second to 10 seconds of CPU, 
and with 11 digits from few tens of seconds to one 
hour. 

If A = L/N is the length of one step, L is the 
length of the whole path and N the total num¬ 
ber of steps, the 4th-order Runge-Kutta formula 
discards terms of order A 5 , so the whole error 
behaves as crk = N A 5 = L 5 /N 4 , and a proper 
choice of L and N allows the control of the pre¬ 
cision. 

Indeed we estimate the relative error, as usual, 
by comparing a value obtained with N steps 
with the one obtained with iV/10 steps, crk = 
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[V(iV) — V'(AT/10)]/V’(AT), to which we add a cu¬ 
mulative rounding error c cre = %/iV x 10” 15 , due 
to our 15 digits double precision FORTRAN im¬ 
plementation. 

The general massive sunrise graph is numer¬ 
ically well studied in literature and several nu¬ 
merical methods are developed, such as multiple 
expansions [12], or numerical integration [12-17]. 
Comparisons are presented in [6] with some val¬ 
ues available in the literature [12,17] with excel¬ 
lent agreement (up to more than 11 digits). 

5. Perspectives 

The presented method for numerically advanc¬ 
ing the solutions of the MDE is rather precise 
and competitive with other available methods for 
numerical MI calculations. 

Rather than conclusions it is more appropriate 
at this stage to present perspectives. It seems 
to be possible to complete the 2-loop self-mass 
for arbitrary internal masses and we have almost 
completed the 4-denominators case [9]. 

We think that the extension to graphs with 
more loops or legs do not present serious prob¬ 
lems, even if the growth in the number of MI in¬ 
creases the computing time. 

It is worth to mention that the method relies on 
the same MDE, which are used also for analytic 
calculations, so it provides a ’low-cost’ comforting 
cross-check for those results. 
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